#include "../rtCommon/F_trace.h"
#include "Random.h"

#include "../rtCommon/Intersect.h"

#include <omp.h>

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

void mpGenerateEyeRays(IN  C_Camera	& camera,
											 IN  int2			& resolution,
											 IN  int2			& sample,
											 IN  uint32			type, 
											 OUT float4		* pRayO,
											 OUT float4		* pRayD )
{
	float3 eye = camera.Eye();
	float3 dir = camera.Dir();
	float3 up  = camera.Up();
	
	float  zn  = camera.zN();
	float  zf  = camera.zF();
	float  asp = camera.Asp();
	float  fov = camera.FoV();

	float3 right = camera.Right();

	// virtual plane
	float2 vp; 
	vp.x = zn * tanf(to_radian(fov * 0.5f));
	vp.y = vp.x / asp;

	float2 scale = { 1.0f, 1.0f };

	if (type == e_Eye_Ray_Perspective)
	{
		scale.x = scale.x / (resolution.x * 0.5f) * vp.x;
		scale.y = scale.y / (resolution.y * 0.5f) * vp.y;

#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (int32 ry = 0; ry < resolution.y; ry++)
			for (int32 rx = 0; rx < resolution.x; rx++)
			{
				float3 Eye;
				float3 Dir;

				for (int32 sy = 0; sy < sample.y; sy++)
					for (int32 sx = 0; sx < sample.x; sx++)
					{
						// sampling plane
						float2 sp = { (sx + 0.5f) / sample.x , (sy + 0.5f) / sample.y };

						float2 pp;
						pp.x = scale.x * (rx - 0.5f * resolution.x + sp.x);
						pp.y = scale.y * (ry - 0.5f * resolution.y + sp.y);

						Dir = normalize(right * pp.x + up * pp.y + dir * zn);

						Eye = eye + Dir;

						uint32 iRay = (sx + rx * sample.x) + (sy + ry * sample.y) * (resolution.x * sample.x);
						pRayD[iRay] = f4(Dir);
						pRayO[iRay] = f4(Eye);
					}
			}	
	}

	if (type == e_Eye_Ray_Orthogonal)
	{
		float3 Dir = dir;

#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (int32 ry = 0; ry < resolution.y; ry++)
			for (int32 rx = 0; rx < resolution.x; rx++)
			{
				for (int32 sy = 0; sy < sample.y; sy++)
					for (int32 sx = 0; sx < sample.x; sx++)
					{
						// sampling plane
						float2 sp = { (sx + 0.5f) / sample.x , (sy + 0.5f) / sample.y };

						float2 pp;
						pp.x = scale.x * (rx - 0.5f * resolution.x + sp.x);
						pp.y = scale.y * (ry - 0.5f * resolution.y + sp.y);

						float3 Eye = eye + right * pp.x + up * pp.y;

						int32 iRay = (sx + rx * sample.x) + (sy + ry * sample.y) * (resolution.x * sample.x);
						pRayD[iRay] = f4(Dir);
						pRayO[iRay] = f4(Eye);
					}
			}
	}
}

void mpGenerateEyeRay(IN  C_Camera & camera,
											IN  int2    & resolution,
											IN  uint32    type, 
											IN  int       x,
											IN  int       y,
											OUT float3  & rayO,
											OUT float3  & rayD,
											OUT int     & iRay)
{
	float3 eye = camera.Eye();
	float3 dir = camera.Dir();
	float3 up  = camera.Up();

	float  zn  = camera.zN();
	float  zf  = camera.zF();
	float  asp = camera.Asp();
	float  fov = camera.FoV();

	float3 right = camera.Right();

	float2 s = { 1.0f, 1.0f };

	// virtual plane
	float2 vp; 
	vp.x = zn * tanf(to_radian(fov * 0.5f));
	vp.y = vp.x / asp;

	int32 ry = y;
	int32 rx = x;

	if (type == e_Eye_Ray_Perspective)
	{
		s.x = s.x / (resolution.x * 0.5f) * vp.x;
		s.y = s.y / (resolution.y * 0.5f) * vp.y;

		float2 sp = {0.5f , 0.5f};

		float2 pp;
		pp.x = s.x * (rx - 0.5f * resolution.x + sp.x);
		pp.y = s.y * (ry - 0.5f * resolution.y + sp.y);
		rayD = right * pp.x + up * pp.y + dir * zn;

		rayD = normalize(rayD);

		rayO = eye + rayD;

		iRay = rx + ry * resolution.x;
	}

	if (type == e_Eye_Ray_Orthogonal)
	{
		rayD = dir;

		float2 sp = {0.5f, 0.5f};

		float2 pp;
		pp.x = s.x * (rx - 0.5f * resolution.x + sp.x);
		pp.y = s.y * (ry - 0.5f * resolution.y + sp.y);

		rayO = eye + right * pp.x + up * pp.y;

		iRay = rx + ry * resolution.x;
	}
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

template<class T>
void mpResolveBuffer_v4(IN  int2     & resolution,
												IN  int2     & sample,
												IN  uint32     filter,	
												OUT T        * pIn,
												OUT T        * pOut)
{
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (int32 ry = 0; ry < resolution.y; ry++)
		for (int32 rx = 0; rx < resolution.x; rx++)
		{
			float4 out = f4();

			for (int32 sy = 0; sy < sample.y; sy++)
				for (int32 sx = 0; sx < sample.x; sx++)
				{
					int32 idx = (sx + rx * sample.x) + (sy + ry * sample.y) * (resolution.x * sample.x);
					out += f4(pIn[idx]);
				}

			int32 idx = rx + ry * (resolution.x);
			out = out * (1.0f / (sample.x * sample.y));

			T result = { out.x, out.y, out.z, out.w };

			pOut[idx] = result;
		}	
}

template<class T>
void mpResolveBuffer_v3(IN  int2     & resolution,
												IN  int2     & sample,
												IN  uint32     filter,	
												OUT T        * pIn,
												OUT T        * pOut)
{
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (int32 ry = 0; ry < resolution.y; ry++)
		for (int32 rx = 0; rx < resolution.x; rx++)
		{
			float3 out = f3();

			for (int32 sy = 0; sy < sample.y; sy++)
				for (int32 sx = 0; sx < sample.x; sx++)
				{
					int32 idx = (sx + rx * sample.x) + (sy + ry * sample.y) * (resolution.x * sample.x);
					out += f3(pIn[idx]);
				}

				int32 idx = rx + ry * (resolution.x);
				out = out * (1.0f / (sample.x * sample.y));

				T result = { out.x, out.y, out.z };

				pOut[idx] = result;
		}	
}

template<class T>
void mpResolveBuffer_v2(IN  int2     & resolution,
												IN  int2     & sample,
												IN  uint32     filter,	
												OUT T        * pIn,
												OUT T        * pOut)
{
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (int32 ry = 0; ry < resolution.y; ry++)
		for (int32 rx = 0; rx < resolution.x; rx++)
		{
			float2 out = f2();

			for (int32 sy = 0; sy < sample.y; sy++)
				for (int32 sx = 0; sx < sample.x; sx++)
				{
					int32 idx = (sx + rx * sample.x) + (sy + ry * sample.y) * (resolution.x * sample.x);
					out += f2(pIn[idx]);
				}

				int32 idx = rx + ry * (resolution.x);
				out = out * (1.0f / (sample.x * sample.y));

				T result = { out.x, out.y };

				pOut[idx] = result;
		}	
}

template<class T>
void mpResolveBuffer_v1(IN  int2     & resolution,
												IN  int2     & sample,
												IN  uint32     filter,	
												OUT T        * pIn,
												OUT T        * pOut)
{
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (int32 ry = 0; ry < resolution.y; ry++)
		for (int32 rx = 0; rx < resolution.x; rx++)
		{
			float1 out = f1();

			for (int32 sy = 0; sy < sample.y; sy++)
				for (int32 sx = 0; sx < sample.x; sx++)
				{
					int32 idx = (sx + rx * sample.x) + (sy + ry * sample.y) * (resolution.x * sample.x);
					out += f1(pIn[idx]);
				}

				int32 idx = rx + ry * (resolution.x);
				out = out * (1.0f / (sample.x * sample.y));

				T result = { out.x };

				pOut[idx] = result;
		}	
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

void mpTraceIndexed(IN  I_SpatialTree * pKdTree,
									  IN  S_Box					& box,
										IN  float4				* pRayO,
										IN  float4				* pRayD, 
										IN  int32						nRay,
										IN  float4				* pVertex,
										IN  int32					* pFace,
										IN  uint32				* pIndex,
										OUT float					* pTmin, 
										OUT int2					* pId, 
										OUT uchar4				* pRgba8)
{
	S_KdNode * pKdNode = (S_KdNode*) (pKdTree->KdTree())->hPtr();

#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (int32 iRay = 0; iRay < nRay; iRay++)
	{
		S_KdNode node_stack[32];
		int node_stack_index = 0;
		int node_index = 0;

		float tnear = 0.0f, tfar = 0.0f;
		float split = 0.0f, tsplit = 0.0f;
		float thit = CU_MAX_F32, tprev = pTmin[iRay];

		int axis = TREE_LEAF;

		float3 rayD = xyz(pRayD[iRay]);
		float3 rayO = xyz(pRayO[iRay]);

		float3 inv_dir = 1.0f / rayD;
		int    sign_dir[3] = { 0, 0, 0 };

		sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
		sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
		sign_dir[2] = (rayD.z >= 0) ? 1 : 0;

		if (Intersect(box.pmin, box.pmax, rayO, rayD, tnear, tfar))
		{
			if (tnear < 0) tnear = 0.0f;
			node_stack[node_stack_index++] = S_KdNode(node_index, tfar);
			tfar = tnear;
		}

		while( node_stack_index > 0 )
		{
			node_stack_index--;
			tnear = tfar;
			tfar = node_stack[node_stack_index].t;
			node_index = node_stack[node_stack_index].next;

			split = pKdNode[node_index].Split();   // as_float(as_int(pNodes[node_index].split) & (~KD_LEAF));
			axis  = pKdNode[node_index].Axis();    // as_int(pNodes[node_index].split) & KD_LEAF;

			tsplit = (split - ((float*)&rayO)[axis]) * ((float*)&inv_dir)[axis];

			int first_left = sign_dir[axis];

			int near = 0;
			int far = 0;

			int left = node_index+1;
			int right = pKdNode[node_index].next;

			if (axis != TREE_LEAF)
			{
				if (first_left)
				{
					near = left;
					far = right;
				}
				else
				{
					far = left;
					near = right;
				}

				if(tsplit < tnear)
				{
					node_index = far;
				}
				else if (tsplit > tfar)
				{
					node_index = near;
				}
				else
				{
					node_index = near;
					node_stack[node_stack_index].t = tfar;
					node_stack[node_stack_index].next = far;

					node_stack_index++;

					tfar = tsplit;
				}
			}
			else
			{
				int32 prim_start = pKdNode[node_index].Start();
				int32 prim_count = pKdNode[node_index].Count();

				for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
				{
					int32 face_idx = pIndex[ prim_start + iPrimitive ];

					int3 face = i3( pFace[3*face_idx], pFace[3*face_idx+1], pFace[3*face_idx+2] );

					float3 v[3] = { xyz(pVertex[face.x]), xyz(pVertex[face.y]), xyz(pVertex[face.z]) };

					if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
					{					
						if (tprev > thit)
						{
							const uchar4 white = {255, 255, 255, 255};

							pRgba8[iRay] = white;
							tprev = thit;					
						} // (tprev > thit)
					}  // if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
				} // for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
			} //if (axis != KD_LEAF)
		} // while (!bhit)
	} // rays
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

void mpTraceIndexed(IN  I_SpatialTree * pKdTree,
									  IN  S_Box					& box,
										IN  float3				& orig,
										IN  float3				& dir,
										IN  int							iRay,
										IN  float4				* pVertex,
										IN  int32					* pIndex,
										IN  uint32				* pTreeIdx,
										OUT float					* pTmin, 
										OUT int2					* pId, 
										OUT uchar4				* pRgba8)
{
	S_KdNode * pKdNode = (S_KdNode *) pKdTree->KdTree();

	S_KdNode node_stack[32];
	int node_stack_index = 0;

	float tnear = 0.0f, tfar = 0.0f;
	float split = 0.0f, tsplit = 0.0f;
	float thit = CU_MAX_F32, tprev = thit;

	int node_index = TREE_ROOT;
	int axis = TREE_LEAF;

	float3 rayO = orig;
	float3 rayD = dir; 

	int2 id = { -1, -1 };

	bool bHit = true;

	float3 inv_dir = 1.0f / rayD;
	int    sign_dir[3] = { 0, 0, 0 };

	sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
	sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
	sign_dir[2] = (rayD.z >= 0) ? 1 : 0;

	if (Intersect(box.pmin, box.pmax, rayO, rayD, tnear, tfar))
	{
		if (tnear < 0) tnear = 0.0f;
		node_stack[node_stack_index++] = S_KdNode(node_index, tfar);
		tfar = tnear;
	}

	while( node_stack_index > 0 )
	{
		node_stack_index--;
		tnear = tfar;
		tfar = node_stack[node_stack_index].t;
		node_index = node_stack[node_stack_index].next;

		split = pKdNode[node_index].Split();
		axis  = pKdNode[node_index].Axis();

		tsplit = (split - ((float*)&rayO)[axis]) * ((float*)&inv_dir)[axis];

		int first_left = sign_dir[axis];

		int near = 0;
		int far = 0;

		int left = node_index+1;
		int right = pKdNode[node_index].next;

		if (axis != TREE_LEAF)
		{
			if (first_left)
			{
				near = left;
				far = right;
			}
			else
			{
				far = left;
				near = right;
			}

			if(tsplit < tnear)
			{
				node_index = far;
			}
			else if (tsplit > tfar)
			{
				node_index = near;
			}
			else
			{
				node_index = near;
				node_stack[node_stack_index].t = tfar;
				node_stack[node_stack_index].next = far;

				node_stack_index++;

				tfar = tsplit;
			}
		}
		else
		{
			int32 prim_start = pKdNode[node_index].Start();
			int32 prim_count = pKdNode[node_index].Count();

			for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
			{
				int32 face_idx = pTreeIdx[ prim_start + iPrimitive ];

				int3 face = i3( pIndex[3*face_idx], pIndex[3*face_idx+1], pIndex[3*face_idx+2] );

				float3 v[3] = { xyz(pVertex[face.x]), xyz(pVertex[face.y]), xyz(pVertex[face.z]) };

				if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
				{					
					if (tprev > thit)
					{
						const uchar4 white = {255, 255, 255, 255};

						pRgba8[iRay] = white;
						tprev = thit;					

					} // (tprev > thit)
				}  // if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
			} // for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
		} //if (axis != KD_LEAF)
	} // while (!bhit)
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

void mpTraceIndexedInstanced(IN  I_SpatialTree * pActorSpatialTree,  
														 IN  I_SpatialTree * pAssetSpatialTree,
														 IN  matrix4x4     * pMatrix,
														 IN  float4       * pRayO, 
														 IN  float4       * pRayD, 
														 IN  int32          nRay,
														 IN  float4       * pPosition, 
														 IN  int32        * pFace, 
														 IN  int32        * pIndex,
														 IN  int32        * pAssetIndex,
														 OUT float        * pTmin, 
														 OUT int2         * pId, 
														 OUT uchar4       * pRgba8)
{
	float4		* pActorData = (float4 * ) pActorSpatialTree->ActorData()->hPtr();
	S_BVHNode * pActorTree = (S_BVHNode *) pActorSpatialTree->BVHTree()->hPtr();
	S_BVHNode * pAssetTree = (S_BVHNode *) pAssetSpatialTree->BVHTree()->hPtr();

#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (int32 iRay = 0; iRay < nRay; iRay++)
	{
		int actor_node_stack[32];
		int actor_node_stack_index = 0;

		float tnear = 0.0f, tfar = 0.0f;
		float split = 0.0f, tsplit = 0.0f;
		float thit = CU_MAX_F32, tprev = pTmin[iRay];

		int node_index = TREE_ROOT;
		int near, far, first_left;

		float3 rayD = xyz(pRayD[iRay]);
		float3 rayO = xyz(pRayO[iRay]);

		float3 inv_dir = 1.0f / rayD;
		int sign_dir[3] = { 0, 0, 0 }; 

		sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
		sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
		sign_dir[2] = (rayD.z >= 0) ? 1 : 0;
		
		actor_node_stack[actor_node_stack_index++] = node_index;

		while (actor_node_stack_index > 0)
		{
			actor_node_stack_index--;

			node_index = actor_node_stack[actor_node_stack_index];
			S_BVHNode bvh = pActorTree[node_index];

			float3 pmin = bvh.pmin;
			float3 pmax = bvh.pmax;

			if (Intersect(pmin, pmax, rayO, rayD, tnear, tfar))
			{
				if (tnear < 0) tnear = 0.0f;

				int axis = bvh.Axis();
				float split = bvh.Split();

				if (axis != TREE_LEAF)
				{
					float o = ((float*)&rayO)[axis];
					float d = ((float*)&inv_dir)[axis];

					tsplit = (split - o) * d;

					first_left = sign_dir[axis];

					if (first_left) 
					{
						near = node_index+1;
						far = bvh.right;
					}
					else
					{
						near = bvh.right;
						far = node_index+1;
					}

					if(tsplit < tnear)
					{
						actor_node_stack[actor_node_stack_index++] = far;
					}
					else if (tsplit > tfar)
					{
						actor_node_stack[actor_node_stack_index++] = near;
					}
					else
					{
						actor_node_stack[actor_node_stack_index++] = far;
						actor_node_stack[actor_node_stack_index++] = near;
					}
				} // if (axis != KD_LEAF)
				else
				{
					int32 actor_start = bvh.Start();
					int32 actor_count = bvh.Count();

					float4 ray_orig_stack[32];
					float4 ray_dir_stack[32];
					int ray_stack_index = 0;

					int asset_node_stack[32];
					int asset_stack_index = 0;

					for (int iActor = 0; iActor < actor_count; iActor++)
					{
						int * pActorTreeIdx = (int *) pActorSpatialTree->TreeIdx()->hPtr();
						int actor_idx = pActorTreeIdx[actor_start+iActor];
						int actor_matrix_idx = pActorData[actor_idx].z;
						int asset_root_idx  = pActorData[actor_idx].y;

						matrix4x4 actor_mat4x4 = pMatrix[actor_matrix_idx];

						Matrix4x4Inverse(&actor_mat4x4, &pMatrix[actor_matrix_idx]);

						ray_orig_stack[ray_stack_index] = f4(rayO);
						ray_dir_stack[ray_stack_index] = f4(rayD);
						ray_stack_index++;

						rayO = Matrix4x4TransformPointR(rayO, actor_mat4x4);
						rayD = normalize(Matrix4x4TransformDirectionR(rayD, actor_mat4x4));

						float3 inv_dir = 1.0f / rayD;
						int sign_dir[3] = { 0, 0, 0 }; 
						sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
						sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
						sign_dir[2] = (rayD.z >= 0) ? 1 : 0;

						asset_node_stack[asset_stack_index++] = asset_root_idx;				

						while (asset_stack_index > 0)
						{
							asset_stack_index--;

							node_index = asset_node_stack[asset_stack_index];
							S_BVHNode bvh = pAssetTree[node_index];

							if (Intersect(bvh.pmin, bvh.pmax, rayO, rayD, tnear, tfar))
							{
								if (tnear < 0) tnear = 0.0f;

								int axis = bvh.Axis();
								float split = bvh.Split();

								if (axis != TREE_LEAF)
								{
									float o = ((float*)&rayO)[axis];
									float d = ((float*)&inv_dir)[axis];

									tsplit = (split - o) * d;

									first_left = sign_dir[axis];

									if (first_left) {
										near = node_index+1;
										far = bvh.right;
									}
									else {
										near = bvh.right;
										far = node_index+1;
									}

									if(tsplit < tnear) asset_node_stack[asset_stack_index++] = far;
									else 
										if (tsplit > tfar) asset_node_stack[asset_stack_index++] = near; 
										else
										{
											asset_node_stack[asset_stack_index++] = far;
											asset_node_stack[asset_stack_index++] = near;
										}
								} // if (axis != KD_LEAF)
								else
								{
									int32 prim_start = bvh.start;
									int32 prim_count = bvh.Count();

									bool bFoundhPrimitive = true;

									for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
									{
										int32 face_idx = pIndex[ prim_start + iPrimitive ];

										int3 face = i3( pFace[3*face_idx], pFace[3*face_idx+1], pFace[3*face_idx+2] );

										float3 v[3] = { xyz(pPosition[face.x]), xyz(pPosition[face.y]), xyz(pPosition[face.z]) };

										if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
										{					
											if (tprev > thit){										
												tprev = thit;					
											  bFoundhPrimitive = false;
											} // if (tprev > thit)
										} //if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
									} // for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)					

									if (!bFoundhPrimitive) {
										const uchar4 white = {255, 255, 255, 255};	
										pRgba8[iRay] = white;
										pTmin[iRay] = tprev;
									} // if (bHit) 
								}
							}
						} // while search for intersection in asset model

						ray_stack_index--;
						rayO = f3(ray_orig_stack[ray_stack_index]);
						rayD = f3(ray_dir_stack[ray_stack_index]);
					}
				} // else (if (axis != KD_LEAF))
			} // if (Intersect(pmin, pmax, rayO, rayD, tnear, tfar))
		} // while (!bHit && stack_index > 0)
	} // rays
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

void mpTraceIndexedInstanced(IN  I_SpatialTree * pActorSpatialTree,  
														 IN  I_SpatialTree * pAssetSpatialTree,
														 IN  matrix4x4     * pMatrix,
														 IN  const float3  & RayO, 
														 IN  const float3  & RayD, 
														 IN  int32					 iRay,
														 IN  float4				* pVertex, 
														 IN  int32        * pIndex, 
														 IN  int32        * pTreeIdx,
														 IN  int32        * pAssetIndex,
														 OUT float        * pTmin, 
														 OUT int2         * pId, 
														 OUT uchar4       * pRgba8)
{
	float4		* pActorData	= (float4 * ) pActorSpatialTree->ActorData()->hPtr();
	S_BVHNode * pActorTree	= (S_BVHNode *) pActorSpatialTree->BVHTree()->hPtr();
	S_BVHNode * pAssetTree	= (S_BVHNode *) pAssetSpatialTree->BVHTree()->hPtr();

	int actor_node_stack[32];
	int actor_node_stack_index = 0;

	float tnear = 0.0f, tfar = 0.0f;
	float split = 0.0f, tsplit = 0.0f;
	float thit = CU_MAX_F32, tprev = pTmin[iRay];

	int node_index = TREE_ROOT;
	int near, far, first_left;

	float3 rayD = RayD;
	float3 rayO = RayO;

	int2 id = { -1, -1 };

	float3 inv_dir = 1.0f / rayD;
	int sign_dir[3] = { 0, 0, 0 }; 

	sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
	sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
	sign_dir[2] = (rayD.z >= 0) ? 1 : 0;

	actor_node_stack[actor_node_stack_index++] = node_index;

	while (actor_node_stack_index > 0)
	{
		actor_node_stack_index--;
		node_index = actor_node_stack[actor_node_stack_index];

		S_BVHNode bvh = pActorTree[node_index];

		if (Intersect(bvh.pmin, bvh.pmax, rayO, rayD, tnear, tfar))
		{
			if (tnear < 0) tnear = 0.0f;

			int axis = bvh.Axis();
			float split = bvh.Split();

			if (axis != TREE_LEAF)
			{
				float o = ((float*)&rayO)[axis];
				float d = ((float*)&inv_dir)[axis];

				tsplit = (split - o) * d;

				first_left = sign_dir[axis];

				if (first_left) 
				{
					near = node_index+1;
					far = pActorTree[node_index].right;
				}
				else
				{
					near = pActorTree[node_index].right;
					far = node_index+1;
				}

				if(tsplit < tnear)
				{
					actor_node_stack[actor_node_stack_index++] = far;
				}
				else if (tsplit > tfar)
				{
					actor_node_stack[actor_node_stack_index++] = near;
				}
				else
				{
					actor_node_stack[actor_node_stack_index++] = far;
					actor_node_stack[actor_node_stack_index++] = near;
				}
			} // if (axis != KD_LEAF)
			else
			{
				int32 actor_start = bvh.Start();
				int32 actor_count = bvh.Count();

				int asset_node_stack[32];
				int asset_stack_index = 0;
				float4 ray_orig_stack[32];
				float4 ray_dir_stack[32];
				int ray_stack_index = 0;

				int * pActorIdx = (int *) pActorSpatialTree->TreeIdx()->hPtr();

				for (int iActor = 0; iActor < actor_count; iActor++)
				{
					int actor_idx = pActorIdx[actor_start+iActor];
					int actor_matrix_idx = pActorData[actor_idx].z;
					int asset_root_idx  = pActorData[actor_idx].y;

					matrix4x4 actor_mat4x4 = pMatrix[actor_matrix_idx];

					Matrix4x4Inverse(&actor_mat4x4, &pMatrix[actor_matrix_idx]);

					ray_orig_stack[ray_stack_index] = f4(rayO);
					ray_dir_stack[ray_stack_index] = f4(rayD);

					rayO = Matrix4x4TransformPointR(rayO, actor_mat4x4);
					rayD = normalize(Matrix4x4TransformDirectionR(rayD, actor_mat4x4));

					float3 inv_dir = 1.0f / rayD;
					int sign_dir[3] = { 0, 0, 0 }; 

					sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
					sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
					sign_dir[2] = (rayD.z >= 0) ? 1 : 0;

					ray_stack_index++;

					asset_node_stack[asset_stack_index++] = asset_root_idx ;

					while (asset_stack_index > 0)
					{
						asset_stack_index--;

						node_index = asset_node_stack[asset_stack_index];
						S_BVHNode bvh = pAssetTree[node_index];

						if (Intersect(bvh.pmin, bvh.pmax, rayO, rayD, tnear, tfar))
						{
							if (tnear < 0) tnear = 0.0f;

							int axis = bvh.Axis();
							float split = bvh.Split();

							if (axis != TREE_LEAF)
							{
								float o = ((float*)&rayO)[axis];
								float d = ((float*)&inv_dir)[axis];

								tsplit = (split - o) * d;

								first_left = sign_dir[axis];

								if (first_left) {
									near = node_index+1;
									far = pAssetTree[node_index].right;
								}
								else {
									near = pAssetTree[node_index].right;
									far = node_index+1;
								}

								if(tsplit < tnear) {
									asset_node_stack[asset_stack_index++] = far;
								}
								else if (tsplit > tfar) {
									asset_node_stack[asset_stack_index++] = near;
								}
								else {
									asset_node_stack[asset_stack_index++] = far;
									asset_node_stack[asset_stack_index++] = near;
								}
							} // if (axis != KD_LEAF)
							else
							{
								int32 prim_start = bvh.start;
								int32 prim_count = bvh.Count();

								bool bFoundPrimitive = false;

								for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
								{
									int32 face_idx = pTreeIdx[ prim_start + iPrimitive ];

									int3 face = i3( pIndex[3*face_idx], pIndex[3*face_idx+1], pIndex[3*face_idx+2] );

									float3 v[3] = { xyz(pVertex[face.x]), xyz(pVertex[face.y]), xyz(pVertex[face.z]) };

									if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
									{					
										if (tprev > thit) {										
											tprev = thit;					
											bFoundPrimitive = true;
										} // if (tprev > thit)
									} //if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
								} // for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)					

								if (bFoundPrimitive) {
									const uchar4 white = {255, 255, 255, 255};	
									pRgba8[iRay] = white;
									pTmin[iRay] = tprev;
								} // if (bHit) 
							}
						}
					} // while search for intersection in asset model

					ray_stack_index--;
					rayO = f3(ray_orig_stack[ray_stack_index]);
					rayD = f3(ray_dir_stack[ray_stack_index]);					
				}
			} // else (if (axis != KD_LEAF))
		} // if (Intersect(pmin, pmax, rayO, rayD, tnear, tfar))
	} // while (!bHit && stack_index > 0)
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

void mpTraceIndexed(IN  I_SpatialTree * pBVHTree,  
										IN  float4				* pRayO, 
										IN  float4				* pRayD,
										IN  int32						nRay,
										IN  float4				* pVertex, 
										IN  int32					* pIndex, 
										IN  uint32				* pTreeIdx,
										OUT float					* pTmin, 
										OUT int2					* pId, 
										OUT uchar4				* pRgba8)
{
	S_BVHNode * pBVHNode = (S_BVHNode *) pBVHTree->BVHTree()->hPtr();

#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (int32 iRay = 0; iRay < nRay; iRay++)
	{
		int node_stack[32];

		float tnear = 0.0f, tfar = 0.0f;
		float split = 0.0f, tsplit = 0.0f;
		float thit = CU_MAX_F32, tprev = pTmin[iRay];
		int stack_index = 0;
		int node_index = TREE_ROOT;
		int near, far, first_left;

		float3 rayD = xyz(pRayD[iRay]);
		float3 rayO = xyz(pRayO[iRay]);

		float3 inv_dir = 1.0f / rayD;
		int sign_dir[3] = { 0, 0, 0 }; 

		sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
		sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
		sign_dir[2] = (rayD.z >= 0) ? 1 : 0;

		node_stack[stack_index++] = node_index;

		while (stack_index > 0)
		{
			stack_index--;
			node_index = node_stack[stack_index];

			S_BVHNode bvh = pBVHNode[node_index];

			if (Intersect(bvh.pmin, bvh.pmax, rayO, rayD, tnear, tfar))
			{
				if (tnear < 0) tnear = 0.0f;

				int axis = bvh.Axis();
				float split = bvh.Split();

				if (axis != TREE_LEAF)
				{
					float o = ((float*)&rayO)[axis];
					float d = ((float*)&inv_dir)[axis];

					tsplit = (split - o) * d;

					first_left = sign_dir[axis];

					if (first_left) {
						near = node_index+1;
						far = pBVHNode[node_index].right;
					}
					else {
						near = pBVHNode[node_index].right;
						far = node_index+1;
					}

					if(tsplit < tnear) {
						node_stack[stack_index++] = far;
					}
					else if (tsplit > tfar)	{
						node_stack[stack_index++] = near;
					}
					else {
						node_stack[stack_index++] = far;
						node_stack[stack_index++] = near;
					}
				} // if (axis != KD_LEAF)
				else
				{
					int32 prim_start = bvh.Start();
					int32 prim_count = bvh.Count();

					for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
					{
						int32 face_idx = pTreeIdx[ prim_start + iPrimitive ];

						int3 face = i3( pIndex[3*face_idx], pIndex[3*face_idx+1], pIndex[3*face_idx+2] );

						float3 v[3] = { xyz(pVertex[face.x]), xyz(pVertex[face.y]), xyz(pVertex[face.z]) };

						if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
						{					
							if (tprev > thit)	{ tprev = thit; } // if (tprev > thit)
						} //if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
					} // for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)					
				} // else (if (axis != KD_LEAF))
			} // if (Intersect(pmin, pmax, rayO, rayD, tnear, tfar))
		} // while (!bHit && stack_index > 0)
	} // rays
}

//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////

void mpTraceIndexed(IN  I_SpatialTree * pBVH,  
										IN  float3				& orig,
										IN  float3				& dir,
										IN  int32						iRay,
										IN  float4				* pVertex, 
										IN  int32					* pIndex, 
										IN  uint32				* pTreeIdx,
										OUT float					* pTmin, 
										OUT int2					* pId, 
										OUT uchar4				* pRgba8)
{
	S_BVHNode * pNodes = (S_BVHNode *) pBVH->BVHTree();

	int node_stack[32];
	int node_stack_index = 0;

	float tnear = 0.0f, tfar = 0.0f;
	float split = 0.0f, tsplit = 0.0f;
	float thit = CU_MAX_F32, tprev = thit;

	int node_index = TREE_ROOT;
	int near, far, first_left;

	float3 rayD = dir;  
	float3 rayO = orig; 

	float3 inv_dir = 1.0f / rayD;
	int sign_dir[3] = { 0, 0, 0 }; 

	sign_dir[0] = (rayD.x >= 0) ? 1 : 0;
	sign_dir[1] = (rayD.y >= 0) ? 1 : 0;
	sign_dir[2] = (rayD.z >= 0) ? 1 : 0;

	node_stack[node_stack_index++] = node_index;

	while (node_stack_index > 0)
	{
		node_stack_index--;
		node_index = node_stack[node_stack_index];

		S_BVHNode bvh = pNodes[node_index];

		if (Intersect(bvh.pmin, bvh.pmax, rayO, rayD, tnear, tfar))
		{
			if (tnear < 0) tnear = 0.0f;

			int axis = bvh.Axis();
			float split = bvh.Split();

			if (axis != TREE_LEAF)
			{
				float o = ((float*)&rayO)[axis];
				float d = ((float*)&inv_dir)[axis];

				tsplit = (split - o) * d;

				first_left = sign_dir[axis];

				if (first_left) 
				{
					near = node_index+1;
					far = pNodes[node_index].right;
				}
				else
				{
					near = pNodes[node_index].right;
					far = node_index+1;
				}

				if(tsplit < tnear) { node_stack[node_stack_index++] = far; }
				else 
					if (tsplit > tfar) { node_stack[node_stack_index++] = near; }
					else
					{
						node_stack[node_stack_index++] = far;
						node_stack[node_stack_index++] = near;
					}
			}
			else
			{
				int32 prim_start = bvh.start;
				int32 prim_count = bvh.Count();

				for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)
				{
					int32 face_idx = pTreeIdx[ prim_start + iPrimitive ];

					int3 face = i3( pIndex[3*face_idx], pIndex[3*face_idx+1], pIndex[3*face_idx+2] );

					float3 v[3] = { xyz(pVertex[face.x]), xyz(pVertex[face.y]), xyz(pVertex[face.z]) };

					if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
					{					
						if (tprev > thit)	{	tprev = thit; } // if (tprev > thit)
					} // if (Intersect(v[0], v[1], v[2], rayO, rayD, thit))
				} // for (int32 iPrimitive = 0; iPrimitive < prim_count; iPrimitive++)					
			} // else (axis != KD_LEAF)
		} // if (Intersect(pmin, pmax, rayO, rayD, tnear, tfar))
	} // while (!bHit && stack_index > 0)
}
